library(Seurat)
library(ggplot2)
library(dplyr)
library(harmony)
library(clusterProfiler)
library(enrichplot)
library(msigdbr)
library(stringr)
library(tidyr)
library(patchwork)
library(Azimuth)
set.seed(2025)

load data

out_path <- "~/projects/2024/RA/Flu_integration/"
# load E_MTAB_12392
in_path <- "~/projects/2024/RA/E_MTAB_12392/out/"
E_MTAB_12392_scdata <- readRDS(paste0(in_path,"E_MTAB_12392_scdata.rds"))
E_MTAB_12392_metadata <- readRDS(paste0(in_path,"E_MTAB_12392_metadata.rds"))
length(E_MTAB_12392_scdata)
[1] 18
# load GSE164381
in_path <- "~/projects/2024/RA/GSE164381_RAW/out/"
GSE164381_scdata <- readRDS(paste0(in_path,"GSE164381_scdata.rds"))
GSE164381_metadata <- readRDS(paste0(in_path,"GSE164381_metadata.rds"))
length(GSE164381_scdata)
[1] 11

select passed cutoff cells in each cohort

# sample name
d_name_1 <- unique(E_MTAB_12392_metadata$orig.ident)
d_name_2 <- unique(GSE164381_metadata$orig.ident)
d_name <- c(d_name_1,d_name_2)
# cell list
cell_list_1 <- rownames(E_MTAB_12392_metadata)
cell_list_2 <- rownames(GSE164381_metadata)
head(cell_list_1)
[1] "S1_AAACCTGAGATAGTCA-1" "S1_AAACCTGAGGAACTGC-1" "S1_AAACCTGGTGTCGCTG-1"
[4] "S1_AAACGGGTCACATGCA-1" "S1_AAACGGGTCCAACCAA-1" "S1_AAACGGGTCCAATGGT-1"
head(cell_list_2)
[1] "GSM5008746_AAACCTGAGTGACATA-1" "GSM5008746_AAACCTGCATTAGCCA-1"
[3] "GSM5008746_AAACCTGGTCTCGTTC-1" "GSM5008746_AAACCTGGTGCCTGGT-1"
[5] "GSM5008746_AAACCTGTCACTTACT-1" "GSM5008746_AAACGGGAGACAGACC-1"
cell_list <- c(cell_list_1, cell_list_2)
# merge
scdata <- c(E_MTAB_12392_scdata, GSE164381_scdata)
Sobj <- merge(scdata[[1]], y = scdata[-1], add.cell.ids = d_name)
Sobj[["percent.mt"]] <- PercentageFeatureSet(Sobj, pattern = "^MT-")
# select passed cells 
Sobj <- subset(Sobj, cells = cell_list)
length(cell_list)
[1] 39055
dim(Sobj)
[1] 36601 39055

normalize, scale, PCA

Sobj <- NormalizeData(Sobj, normalization.method = "LogNormalize", scale.factor = 10000)
Normalizing layer: counts.S1
Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Normalizing layer: counts.S3
Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Normalizing layer: counts.S11
Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Normalizing layer: counts.S13
Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Normalizing layer: counts.S14
Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Normalizing layer: counts.S16
Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Normalizing layer: counts.S17
Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Normalizing layer: counts.S18
Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Normalizing layer: counts.S19
Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Normalizing layer: counts.S20
Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Normalizing layer: counts.S21
Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Normalizing layer: counts.S22
Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Normalizing layer: counts.S23
Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Normalizing layer: counts.S24
Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Normalizing layer: counts.S25
Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Normalizing layer: counts.S26
Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Normalizing layer: counts.S27
Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Normalizing layer: counts.S28
Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Normalizing layer: counts.GSM5008746
Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Normalizing layer: counts.GSM5008747
Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Normalizing layer: counts.GSM5008748
Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Normalizing layer: counts.GSM5008749
Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Normalizing layer: counts.GSM5008750
Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Normalizing layer: counts.GSM5008751
Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Normalizing layer: counts.GSM5008752
Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Normalizing layer: counts.GSM5008753
Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Normalizing layer: counts.GSM5008754
Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Normalizing layer: counts.GSM5008755
Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Normalizing layer: counts.GSM5008756
Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Sobj <- FindVariableFeatures(Sobj, selection.method = "vst", nfatures = 2000)
Finding variable features for layer counts.S1
Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Finding variable features for layer counts.S3
Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Finding variable features for layer counts.S11
Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Finding variable features for layer counts.S13
Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  at  -0.31854
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  radius  0.00030656
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  all data on boundary of neighborhood. make span bigger
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  pseudoinverse used at -0.31854
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  neighborhood radius 0.017509
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  reciprocal condition number  1
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  zero-width neighborhood. make span bigger
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  There are other near singularities as well. 0.090619
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Finding variable features for layer counts.S14
Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Finding variable features for layer counts.S16
Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Finding variable features for layer counts.S17
Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Finding variable features for layer counts.S18
Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Finding variable features for layer counts.S19
Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Finding variable features for layer counts.S20
Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  pseudoinverse used at -1.9345
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  neighborhood radius 0.30103
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  reciprocal condition number  6.7402e-15
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Finding variable features for layer counts.S21
Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  pseudoinverse used at -1.752
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  neighborhood radius 0.30103
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  reciprocal condition number  1.0124e-14
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Finding variable features for layer counts.S22
Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  at  -0.91849
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  radius  0.00023715
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  all data on boundary of neighborhood. make span bigger
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  pseudoinverse used at -0.91849
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  neighborhood radius 0.0154
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  reciprocal condition number  1
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  zero-width neighborhood. make span bigger
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  There are other near singularities as well. 0.031008
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Finding variable features for layer counts.S23
Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Finding variable features for layer counts.S24
Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  pseudoinverse used at -2.1206
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  neighborhood radius 0.30103
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  reciprocal condition number  8.0208e-15
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Finding variable features for layer counts.S25
Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  pseudoinverse used at -1.2225
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  neighborhood radius 0.31945
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  reciprocal condition number  7.6188e-28
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  There are other near singularities as well. 0.031008
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Finding variable features for layer counts.S26
Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  at  -0.62058
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  radius  0.00034289
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  all data on boundary of neighborhood. make span bigger
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  pseudoinverse used at -0.62058
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  neighborhood radius 0.018517
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  reciprocal condition number  1
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  zero-width neighborhood. make span bigger
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  There are other near singularities as well. 0.031008
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Finding variable features for layer counts.S27
Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  at  -1.0144
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  radius  0.0002087
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  all data on boundary of neighborhood. make span bigger
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  pseudoinverse used at -1.0144
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  neighborhood radius 0.014447
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  reciprocal condition number  1
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  zero-width neighborhood. make span bigger
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  There are other near singularities as well. 0.031008
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Finding variable features for layer counts.S28
Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  pseudoinverse used at -2.2407
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  neighborhood radius 0.49766
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  reciprocal condition number  1.8412e-14
Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,  :
  There are other near singularities as well. 0.090619
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Finding variable features for layer counts.GSM5008746
Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Finding variable features for layer counts.GSM5008747
Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Finding variable features for layer counts.GSM5008748
Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Finding variable features for layer counts.GSM5008749
Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Finding variable features for layer counts.GSM5008750
Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Finding variable features for layer counts.GSM5008751
Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Finding variable features for layer counts.GSM5008752
Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Finding variable features for layer counts.GSM5008753
Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Finding variable features for layer counts.GSM5008754
Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Finding variable features for layer counts.GSM5008755
Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Finding variable features for layer counts.GSM5008756
Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Sobj <- ScaleData(Sobj)
Centering and scaling data matrix

  |                                                                                             
  |                                                                                       |   0%
  |                                                                                             
  |============================================                                           |  50%
  |                                                                                             
  |=======================================================================================| 100%
Sobj <- RunPCA(Sobj)
PC_ 1 
Positive:  MT-CO2, MT-CO3, MT-CYB, CD83, MALAT1, MS4A1, MT-CO1, CD74, MT-ATP6, SNX9 
       BACH2, HLA-DQB1, ZNF331, SLC2A3, CD79A, KLF6, HLA-DRA, FCER2, MT-ND5, CD69 
       PHACTR1, HLA-DRB1, LY9, IGHM, NEIL1, IL4R, SERPINB9, CCR7, HLA-DQA1, IGHD 
Negative:  GAPDH, ACTG1, S100A4, PFN1, SRGN, ANXA1, ACTB, S100A6, AIF1, CLIC1 
       HCST, TPI1, CTSW, SH3BGRL3, TXN, TALDO1, TYROBP, CFL1, GSTP1, IFITM2 
       ENO1, IGFBP7, FCER1G, VIM, S100A10, MYL12A, PRSS57, CD99, PPIB, ARPC1B 
PC_ 2 
Positive:  HLA-DRA, CD74, HLA-DPA1, HLA-DPB1, RPS8, RPLP0, HLA-DRB1, RPLP1, EEF1B2, RACK1 
       RPS18, RPS4X, RPL5, RPS12, RPS6, RPS2, EEF1A1, RPL7A, RPL10, RPS3A 
       CD79A, RPS14, RPS24, FTL, RPS3, HLA-DQB1, SNHG29, RPL7, HLA-DRB5, FTH1 
Negative:  GNLY, NKG7, CST7, GZMB, CCL5, GZMA, CTSW, CD247, CD7, GZMM 
       FGFBP2, PRF1, KLRD1, GZMH, SPON2, KLRB1, TYROBP, CLIC3, CCL4, HCST 
       ZAP70, FCGR3A, IL2RB, IFITM1, IL32, S1PR5, KLRF1, TRBC1, PRKCH, MATK 
PC_ 3 
Positive:  PRSS57, H1F0, CYTL1, SMIM24, EGFL7, GATA2, AL157895.1, BEX3, GIHCG, NFE2 
       SLC40A1, CPA3, SOX4, MT-CO1, MT-CO3, CD34, EREG, EMID1, MT-CYB, FCER1A 
       MPP1, CRHBP, LAPTM4B, CDK6, IL18, ERG, DEPTOR, LINC02573, NPR3, MALAT1 
Negative:  TMSB4X, HLA-DPA1, TMSB10, B2M, CD79A, ACTB, RPS18, LTB, FTL, RPL10 
       HLA-DPB1, CYBA, EEF1A1, RPLP1, RPS3, CD74, RPS12, HLA-DRA, LSP1, CORO1A 
       RPL7A, HLA-DRB1, CFL1, JUNB, RPS4X, PFN1, RPS14, RPS6, EEF1B2, DUSP1 
PC_ 4 
Positive:  RPS3, EEF1A1, RPS12, RPL10, RPL7A, RPS3A, RPS4X, RPS14, RPL5, RPS8 
       RACK1, RPS6, RPS18, EEF1B2, RPS2, RPLP0, TPT1, RPLP1, RPS24, CTSW 
       NKG7, RPL7, CST7, NPM1, SNHG29, GAS5, CD7, GNLY, HNRNPA1, GZMA 
Negative:  LYZ, CST3, IFI30, S100A9, CLEC7A, FCN1, CPVL, SERPINA1, S100A8, CSTA 
       PLAUR, LGALS3, CFP, LGALS2, CFD, S100A11, TNFAIP2, CD68, G0S2, MAFB 
       AIF1, S100A12, TGFBI, MS4A6A, VCAN, TYMP, SHTN1, AC020656.1, HMOX1, PLBD1 
PC_ 5 
Positive:  TXNDC5, JCHAIN, TNFRSF17, MZB1, IGHA1, CD27, AQP3, PDIA4, DERL3, IGHG1 
       IGHA2, ITM2C, SEC11C, MTRNR2L12, FKBP11, IGKC, IGHG2, PRDM1, SRM, CHPF 
       CLPTM1L, MAN1A1, SDF2L1, CD38, MANF, LMAN1, PRDX4, STT3A, SSR3, GPRC5D 
Negative:  FTH1, LYZ, CST3, RPS14, RPS3A, RPL10, S100A9, FCN1, RPS12, RPS8 
       RPS2, EEF1A1, TYROBP, AIF1, TPT1, RPS24, FCER1G, SERPINA1, S100A8, CSTA 
       CLEC7A, FTL, RPLP1, RPL7A, RPS4X, RPS3, RPL5, PLAUR, LST1, CPVL 

select cutoff of dimension.

ElbowPlot(Sobj)

PCA space

DimPlot(Sobj, group.by = "orig.ident")

UMAP

Sobj <- RunUMAP(Sobj, dims = 1:15, reduction = "pca", reduction.name = "umap.unintegrated")
20:06:19 UMAP embedding parameters a = 0.9922 b = 1.112
20:06:19 Read 39055 rows and found 15 numeric columns
20:06:19 Using Annoy for neighbor search, n_neighbors = 30
20:06:19 Building Annoy index with metric = cosine, n_trees = 50
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
20:06:24 Writing NN index file to temp file /tmp/RtmpBxlMZd/file5c7036bb70824
20:06:24 Searching Annoy index using 1 thread, search_k = 3000
20:06:38 Annoy recall = 100%
20:06:39 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
20:06:43 Initializing from normalized Laplacian + noise (using RSpectra)
20:06:46 Commencing optimization for 200 epochs, with 1667658 positive edges
20:06:46 Using rng type: pcg
Using method 'umap'
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
20:07:08 Optimization finished
p0 <- DimPlot(Sobj, reduction = "umap.unintegrated", group.by = c("orig.ident", "Batch"))
Warning: The following requested variables were not found: Batch
p0

Harmony batch correction

Sobj <- RunHarmony(Sobj,"orig.ident", plot_convergence = T, theta = 3)
Transposing data matrix
Initializing state using k-means centroids initialization
Harmony 1/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 2/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 3/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 4/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony 5/10
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Harmony converged after 5 iterations

harmony_embeddings <- Embeddings(Sobj, "harmony")
harmony_embeddings[1:5, 1:5]
                       harmony_1 harmony_2  harmony_3   harmony_4   harmony_5
S1_AAACCTGAGATAGTCA-1   4.701568 -1.708613  1.3231752  -2.8416889  -0.1320850
S1_AAACCTGAGGAACTGC-1   2.451071  1.390329 -1.8605446   1.6223532   0.2989081
S1_AAACCTGGTGTCGCTG-1   4.284953  2.023095 -0.8945279   1.5819618  -1.1558085
S1_AAACGGGTCACATGCA-1   3.936914  1.105024 -1.2764343   0.8342876  -1.0429716
S1_AAACGGGTCCAACCAA-1 -23.461463  1.001002 -3.2267307 -25.5267511 -14.3578007
DimPlot(Sobj, reduction = "harmony", group.by = "orig.ident")

UMAP with harmony

Sobj <- RunUMAP(Sobj, dims = 1:15, reduction = "harmony")
20:07:59 UMAP embedding parameters a = 0.9922 b = 1.112
20:07:59 Read 39055 rows and found 15 numeric columns
20:07:59 Using Annoy for neighbor search, n_neighbors = 30
20:07:59 Building Annoy index with metric = cosine, n_trees = 50
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
20:08:04 Writing NN index file to temp file /tmp/RtmpBxlMZd/file5c7034015e3f0
20:08:04 Searching Annoy index using 1 thread, search_k = 3000
20:08:19 Annoy recall = 100%
20:08:20 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
20:08:23 Initializing from normalized Laplacian + noise (using RSpectra)
20:08:26 Commencing optimization for 200 epochs, with 1710788 positive edges
20:08:26 Using rng type: pcg
Using method 'umap'
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
20:08:48 Optimization finished
DimPlot(Sobj, reduction = "umap", group.by = "orig.ident")

Check whether the harmony batch correction work well or not.

Idents(Sobj) <- Sobj$orig.ident
rename_map <- c(
    setNames(rep("E_MTAB_12392", length(d_name_1)), d_name_1),
    setNames(rep("GSE164381", length(d_name_2)), d_name_2)
)
Sobj <- RenameIdents(Sobj, rename_map)
unique(Idents(Sobj))
[1] E_MTAB_12392 GSE164381   
Levels: E_MTAB_12392 GSE164381
Sobj$batch <- Idents(Sobj)
DimPlot(Sobj, group.by = "batch")

DimPlot(Sobj, reduction = "umap", group.by = "batch")

Sobj
An object of class Seurat 
36607 features across 39055 samples within 2 assays 
Active assay: RNA (36601 features, 2000 variable features)
 59 layers present: counts.S1, counts.S3, counts.S11, counts.S13, counts.S14, counts.S16, counts.S17, counts.S18, counts.S19, counts.S20, counts.S21, counts.S22, counts.S23, counts.S24, counts.S25, counts.S26, counts.S27, counts.S28, counts.GSM5008746, counts.GSM5008747, counts.GSM5008748, counts.GSM5008749, counts.GSM5008750, counts.GSM5008751, counts.GSM5008752, counts.GSM5008753, counts.GSM5008754, counts.GSM5008755, counts.GSM5008756, data.S1, data.S3, data.S11, data.S13, data.S14, data.S16, data.S17, data.S18, data.S19, data.S20, data.S21, data.S22, data.S23, data.S24, data.S25, data.S26, data.S27, data.S28, data.GSM5008746, data.GSM5008747, data.GSM5008748, data.GSM5008749, data.GSM5008750, data.GSM5008751, data.GSM5008752, data.GSM5008753, data.GSM5008754, data.GSM5008755, data.GSM5008756, scale.data
 1 other assay present: HTO
 4 dimensional reductions calculated: pca, umap.unintegrated, harmony, umap
Sobj <- FindNeighbors(Sobj, dims = 1:15, reduction = "harmony")
Computing nearest neighbor graph
Computing SNN
Sobj <- FindClusters(Sobj, reduction = "umap", resolution = 1)
Warning: The following arguments are not used: reduction
Warning: The following arguments are not used: reduction
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 39055
Number of edges: 1156538

Running Louvain algorithm...
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Maximum modularity in 10 random starts: 0.8409
Number of communities: 18
Elapsed time: 10 seconds
Sobj@reductions
$pca
A dimensional reduction object with key PC_ 
 Number of dimensions: 50 
 Number of cells: 39055 
 Projected dimensional reduction calculated:  FALSE 
 Jackstraw run: FALSE 
 Computed using assay: RNA 

$umap.unintegrated
A dimensional reduction object with key umapunintegrated_ 
 Number of dimensions: 2 
 Number of cells: 39055 
 Projected dimensional reduction calculated:  FALSE 
 Jackstraw run: FALSE 
 Computed using assay: RNA 

$harmony
A dimensional reduction object with key harmony_ 
 Number of dimensions: 50 
 Number of cells: 39055 
 Projected dimensional reduction calculated:  TRUE 
 Jackstraw run: FALSE 
 Computed using assay: RNA 

$umap
A dimensional reduction object with key umap_ 
 Number of dimensions: 2 
 Number of cells: 39055 
 Projected dimensional reduction calculated:  FALSE 
 Jackstraw run: FALSE 
 Computed using assay: RNA 
DimPlot(Sobj, reduction = "umap", group.by = "RNA_snn_res.1")

Sobj <- JoinLayers(Sobj)
DimPlot(Sobj, reduction = "umap", label = T)

load reference data

SLE_Sobj <- readRDS("~/projects/2024/RA/Integration_SLE_Sobj.rds")
SLE_Sobj <- RunUMAP(SLE_Sobj, reduction = "harmony",dims = 1:15, return.model = TRUE)
UMAP will return its model
20:10:08 UMAP embedding parameters a = 0.9922 b = 1.112
20:10:08 Read 133389 rows and found 15 numeric columns
20:10:08 Using Annoy for neighbor search, n_neighbors = 30
20:10:08 Building Annoy index with metric = cosine, n_trees = 50
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
20:10:24 Writing NN index file to temp file /tmp/RtmpBxlMZd/file5c703562f00f7
20:10:24 Searching Annoy index using 1 thread, search_k = 3000
20:11:27 Annoy recall = 100%
20:11:29 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
20:11:34 Initializing from normalized Laplacian + noise (using RSpectra)
20:11:47 Commencing optimization for 200 epochs, with 5832214 positive edges
20:11:47 Using rng type: pcg
Using method 'umap'
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
20:13:09 Optimization finished
DimPlot(SLE_Sobj, reduction = "umap", label = T, repel = T)
Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

DimPlot(SLE_Sobj, reduction = "umap.unintegrated")
Rasterizing points since number of points exceeds 100,000.
To disable this behavior set `raster=FALSE`

# if (!file.exists(paste0(out_path,"Integration_Virus.RData"))) {
#   save.image(file = paste0(out_path,"Integration_Virus.RData"))
# }

calculate transfer anchors

anchors <- FindTransferAnchors(reference = SLE_Sobj,
                               query = Sobj,
                               normalization.method = "LogNormalize",
                               reference.reduction = "harmony",
                               dims = 1:15)
Projecting cell embeddings
Finding neighborhoods
Finding anchors
    Found 6697 anchors

Mapping

Sobj <- MapQuery(anchorset = anchors,
                 query = Sobj,
                 reference = SLE_Sobj,
                 refdata = list(celltype = "anno_text"),
                 reference.reduction = "harmony",
                 reduction.model = "umap")
Finding integration vectors
Finding integration vector weights
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Predicting cell labels
Warning: Layer counts isn't present in the assay object; returning NULL
Warning: Layer counts isn't present in the assay object; returning NULL
Warning: Layer counts isn't present in the assay object; returning NULL

  |                                                  | 0 % ~calculating  

Integrating dataset 2 with reference dataset
Finding integration vectors
Integrating data

  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=02s  
Computing nearest neighbors
Running UMAP projection
20:16:20 Read 39055 rows
20:16:20 Processing block 1 of 1
20:16:20 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
20:16:21 Initializing by weighted average of neighbor coordinates using 1 thread
20:16:21 Commencing optimization for 67 epochs, with 1171650 positive edges
Using method 'umap'
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
20:16:30 Finished
Idents(Sobj) <- Sobj$predicted.celltype
Idents(Sobj) <- factor(x = Idents(Sobj), levels = c("Naive B", "Activated B", "Unswitched memory B", "CD27+ memory B", "Plasmablast", "Plasma", "T cell like B", "Myeloid like B", "Megakaryocyte like B"))
levels(Idents(Sobj))
[1] "Naive B"              "Activated B"          "Unswitched memory B"  "CD27+ memory B"      
[5] "Plasmablast"          "Plasma"               "T cell like B"        "Myeloid like B"      
[9] "Megakaryocyte like B"
levels(SLE_Sobj$anno_text)
 [1] "Naive B"              "Activated B"          "Unswitched memory B"  "CD27+ memory B"      
 [5] "Atypical memory B"    "INF-induced memory B" "Plasmablast"          "Plasma"              
 [9] "T cell like B"        "Myeloid like B"       "Megakaryocyte like B"
Sobj$predicted.celltype <- Idents(Sobj)
table(Sobj$predicted.celltype)

             Naive B          Activated B  Unswitched memory B       CD27+ memory B 
               29127                 2631                   53                 1264 
         Plasmablast               Plasma        T cell like B       Myeloid like B 
                  21                  387                 3885                 1684 
Megakaryocyte like B 
                   3 
Idents(Sobj) <- Sobj$predicted.celltype
Sobj <- subset(Sobj, idents = setdiff(levels(Sobj), "Megakaryocyte like B"))
table(Sobj$predicted.celltype)

            Naive B         Activated B Unswitched memory B      CD27+ memory B 
              29127                2631                  53                1264 
        Plasmablast              Plasma       T cell like B      Myeloid like B 
                 21                 387                3885                1684 
Idents(Sobj) <- Sobj$RNA_snn_res.1
# p1 <- DimPlot(Sobj, reduction = "umap", group.by = "RNA_snn_res.1")
p2 <- DimPlot(Sobj, reduction = "umap", label = T)
Idents(Sobj) <- Sobj$predicted.celltype
p3 <- DimPlot(Sobj, label = TRUE, repel = TRUE, reduction = "ref.umap") + NoLegend()
p4 <- DimPlot(Sobj, label = TRUE, repel = TRUE, reduction = "ref.umap")
p5 <- DimPlot(Sobj, label = TRUE, repel = TRUE, reduction = "umap") + NoLegend()
p6 <- DimPlot(Sobj, label = TRUE, repel = TRUE, reduction = "umap")
E_MTAB_12392_metadata$anno_ABC_LRB <- NA

E_MTAB_12392_metadata$anno_ABC_LRB[E_MTAB_12392_metadata$LSP1_score_pos==T & E_MTAB_12392_metadata$ABC_score_pos==F] <- "LRB"
E_MTAB_12392_metadata$anno_ABC_LRB[E_MTAB_12392_metadata$LSP1_score_pos==F & E_MTAB_12392_metadata$ABC_score_pos==T] <- "ABC"
E_MTAB_12392_metadata$anno_ABC_LRB[E_MTAB_12392_metadata$LSP1_score_pos==T & E_MTAB_12392_metadata$ABC_score_pos==T] <- "LRB & ABC"
E_MTAB_12392_metadata$anno_ABC_LRB[E_MTAB_12392_metadata$LSP1_score_pos==F & E_MTAB_12392_metadata$ABC_score_pos==F] <- "Classical B"

unique(E_MTAB_12392_metadata$anno_ABC_LRB)
[1] "Classical B" "ABC"         "LRB & ABC"   "LRB"        
unique(GSE164381_metadata$anno_ABC_LRB)
[1] ABC         Classical B LRB         LRB & ABC  
Levels: LRB ABC LRB & ABC Classical B
anno_vars <- c("anno_ABC_LRB", "LSP1_score", "ABC_score", "LSP1_score_pos", "ABC_score_pos")

for (var in anno_vars) {
  l1 <- E_MTAB_12392_metadata[[var]]
  l2 <- GSE164381_metadata[[var]]
  
  if (is.factor(l1)) l1 <- as.character(l1)
  if (is.factor(l2)) l2 <- as.character(l2)
  
  names(l1) <- rownames(E_MTAB_12392_metadata)
  names(l2) <- rownames(GSE164381_metadata)
  
  combined_anno <- c(l1, l2)
  
  Sobj[[var]] <- combined_anno[colnames(Sobj)]
}
rm(l1, l2, combined_anno)
my_colors <- c("LRB" = "#ff0000", "ABC" = "#ffff00", "LRB & ABC" = "#00b0ff", "Classical B" = "#BFBFBF")
Idents(Sobj) <- Sobj$anno_ABC_LRB
p7 <- DimPlot(Sobj, reduction = "umap", cols = my_colors) + NoLegend()
DimPlot(Sobj, reduction = "umap", cols = my_colors)

celltype.prop <- as.data.frame(table(Sobj$predicted.celltype, Sobj$anno_ABC_LRB))
# head(celltype.prop)
celltype.prop <- celltype.prop %>%
    rename(celltype = Var1, `ABC LRB` = Var2) %>%
    group_by(celltype) %>%
    mutate(Percent = Freq / sum(Freq)*100)
ggplot(celltype.prop, aes(x = celltype, y = Percent, fill = `ABC LRB`))+
    geom_bar(stat = "identity", width = 0.9) + 
    scale_fill_manual(values = my_colors) +
    theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1), axis.title.x = element_blank(), plot.margin = margin(10, 10, 20, 10))



celltype_percent_wide <- celltype.prop %>%
  select(celltype, `ABC LRB`, Percent) %>%
  pivot_wider(names_from = `ABC LRB`, values_from = Percent)
head(celltype_percent_wide)

celltype_freq_wide <- celltype.prop %>%
  select(celltype, `ABC LRB`, Freq) %>%
  pivot_wider(names_from = `ABC LRB`, values_from = Freq)
head(celltype_freq_wide)
write.csv(celltype_percent_wide, paste0(out_path,"out/virus_barplot_percent.csv"))
write.csv(celltype_freq_wide, paste0(out_path,"out/virus_barplot_freq.csv"))
p2

p3

p4

p5

p6

p7
pdf(paste0(out_path, "Azimuth_predicted_umap_063025.pdf"), width = 4, height = 4)
p2
p3
p4
p5
p6
p7
dev.off()
png 
  2 

Idents(Sobj) <- Sobj$predicted.celltype
p6 <- DimPlot(Sobj, label = TRUE, repel = TRUE, reduction = "umap")
pdf(paste0(out_path, "new_virus_umap_080425.pdf"), width = 4, height = 4)
p6
p7
dev.off()
null device 
          1 
p6_1 <- DimPlot(Sobj, label = FALSE, repel = TRUE, reduction = "umap") +
  labs(x = "", y = "") +
  theme(
    legend.position = "none",
    axis.text.x = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks.x = element_line(),
    axis.ticks.y = element_line()
    )
tiff(filename = paste0(out_path,"umap_virus.tif"), width = 4, height = 4, units = "in", res = 300, compression = "lzw")
print(p6_1)
dev.off()
null device 
          1 
p7_1 <- p7 +
  labs(x = "", y = "") +
  theme(
    legend.position = "none",
    axis.text.x = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks.x = element_line(),
    axis.ticks.y = element_line()
    )
tiff(filename = paste0(out_path,"umap_LRB_ABC.tif"), width = 4, height = 4, units = "in", res = 300, compression = "lzw")
print(p7_1)
dev.off()
null device 
          1 
table(Sobj$anno_ABC_LRB)

        ABC Classical B         LRB   LRB & ABC 
       8486       18588        6126        5852 
library(ggpointdensity)
library(scales)
df <- FetchData(Sobj, vars = c("umap_1", "umap_2", "LSP1_score_pos", "ABC_score_pos", "anno_ABC_LRB"))
df_LSP1 <- df %>% filter(LSP1_score_pos)
df_ABC <- df %>% filter(ABC_score_pos)
t1 <- ggplot(df_LSP1, aes(x = umap_1, y = umap_2)) +
  geom_pointdensity(method = "neighbors", size = 0.1) +
  scale_color_gradientn(
    colors = c("grey", "grey", "red"),
    values = rescale(c(0, 350, 1200)),
    limits = c(0, 1000)
  ) +
    theme_classic()

t2 <- ggplot(df_ABC, aes(x = umap_1, y = umap_2)) +
  geom_pointdensity(method = "neighbors", size = 0.1) +
  scale_color_gradientn(
    colors = c("grey", "grey", "red"),
    values = rescale(c(0, 350, 1200)),
    limits = c(0, 1000)
  ) +
    theme_classic()
t1

t2

pdf(paste0(out_path,"LRB_ABC_umap_density_080425.pdf"), width = 6, height = 5)
t1
t2
dev.off()
null device 
          1 
t1_1 <- t1 +
labs(x = "", y = "") +
theme(
  legend.position = "none",
  axis.text.x = element_blank(),
  axis.text.y = element_blank(),
  axis.ticks.x = element_line(),
  axis.ticks.y = element_line()
  )
t2_1 <- t2 +
labs(x = "", y = "") +
theme(
  legend.position = "none",
  axis.text.x = element_blank(),
  axis.text.y = element_blank(),
  axis.ticks.x = element_line(),
  axis.ticks.y = element_line()
  )
t1_1

t2_1

tiff(filename = paste0(out_path,"LRB_density_umap.tif"), width = 4, height = 4, units = "in", res = 300, compression = "lzw")
print(t1_1)
dev.off()
null device 
          1 
tiff(filename = paste0(out_path,"ABC_density_umap.tif"), width = 4, height = 4, units = "in", res = 300, compression = "lzw")
print(t2_1)
dev.off()
null device 
          1 
if (!file.exists(paste0(out_path,"virus_integrated_data_080425.rds"))) {
  saveRDS(Sobj, file = paste0(out_path,"virus_integrated_data_080425.rds"))
}
sessionInfo()
R version 4.4.2 (2024-10-31)
Platform: x86_64-pc-linux-gnu
Running under: Ubuntu 20.04.2 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0 
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=en_US.UTF-8       
 [4] LC_COLLATE=en_US.UTF-8     LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                  LC_ADDRESS=C              
[10] LC_TELEPHONE=C             LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

time zone: America/Los_Angeles
tzcode source: system (glibc)

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] scales_1.4.0           ggpointdensity_0.2.0   future_1.58.0          Azimuth_0.5.0         
 [5] shinyBS_0.61.1         patchwork_1.3.1        tidyr_1.3.1            stringr_1.5.1         
 [9] msigdbr_25.1.0         enrichplot_1.26.6      clusterProfiler_4.14.6 harmony_1.2.3         
[13] Rcpp_1.1.0             dplyr_1.1.4            ggplot2_3.5.2          Seurat_5.3.0          
[17] SeuratObject_5.1.0     sp_2.2-0              

loaded via a namespace (and not attached):
  [1] fs_1.6.6                          ProtGenerics_1.38.0              
  [3] matrixStats_1.5.0                 spatstat.sparse_3.1-0            
  [5] bitops_1.0-9                      DirichletMultinomial_1.48.0      
  [7] TFBSTools_1.44.0                  httr_1.4.7                       
  [9] RColorBrewer_1.1-3                tools_4.4.2                      
 [11] sctransform_0.4.2                 R6_2.6.1                         
 [13] DT_0.33                           lazyeval_0.2.2                   
 [15] uwot_0.2.3                        rhdf5filters_1.18.1              
 [17] withr_3.0.2                       gridExtra_2.3                    
 [19] progressr_0.15.1                  cli_3.6.5                        
 [21] Biobase_2.66.0                    spatstat.explore_3.4-3           
 [23] fastDummies_1.7.5                 EnsDb.Hsapiens.v86_2.99.0        
 [25] shinyjs_2.1.0                     labeling_0.4.3                   
 [27] sass_0.4.10                       spatstat.data_3.1-6              
 [29] readr_2.1.5                       ggridges_0.5.6                   
 [31] pbapply_1.7-2                     Rsamtools_2.22.0                 
 [33] yulab.utils_0.2.0                 gson_0.1.0                       
 [35] DOSE_4.0.1                        R.utils_2.13.0                   
 [37] parallelly_1.45.0                 BSgenome_1.74.0                  
 [39] rstudioapi_0.14                   RSQLite_2.4.1                    
 [41] generics_0.1.4                    gridGraphics_0.5-1               
 [43] BiocIO_1.16.0                     gtools_3.9.5                     
 [45] ica_1.0-3                         spatstat.random_3.4-1            
 [47] googlesheets4_1.1.1               GO.db_3.20.0                     
 [49] Matrix_1.7-3                      S4Vectors_0.44.0                 
 [51] abind_1.4-8                       R.methodsS3_1.8.2                
 [53] lifecycle_1.0.4                   yaml_2.3.10                      
 [55] SummarizedExperiment_1.36.0       rhdf5_2.50.2                     
 [57] qvalue_2.38.0                     SparseArray_1.6.2                
 [59] Rtsne_0.17                        grid_4.4.2                       
 [61] blob_1.2.4                        promises_1.3.3                   
 [63] shinydashboard_0.7.3              pwalign_1.2.0                    
 [65] crayon_1.5.3                      miniUI_0.1.2                     
 [67] ggtangle_0.0.7                    lattice_0.20-45                  
 [69] cowplot_1.2.0                     annotate_1.84.0                  
 [71] GenomicFeatures_1.58.0            KEGGREST_1.46.0                  
 [73] pillar_1.11.0                     knitr_1.50                       
 [75] fgsea_1.32.4                      GenomicRanges_1.58.0             
 [77] rjson_0.2.23                      future.apply_1.20.0              
 [79] codetools_0.2-19                  fastmatch_1.1-6                  
 [81] glue_1.8.0                        ggfun_0.1.9                      
 [83] spatstat.univar_3.1-3             data.table_1.17.6                
 [85] vctrs_0.6.5                       png_0.1-8                        
 [87] treeio_1.30.0                     spam_2.11-1                      
 [89] cellranger_1.1.0                  poweRlaw_1.0.0                   
 [91] gtable_0.3.6                      assertthat_0.2.1                 
 [93] cachem_1.1.0                      xfun_0.52                        
 [95] Signac_1.14.0                     S4Arrays_1.6.0                   
 [97] mime_0.13                         survival_3.5-0                   
 [99] gargle_1.5.2                      RcppRoll_0.3.1                   
[101] fitdistrplus_1.2-4                ROCR_1.0-11                      
[103] nlme_3.1-162                      ggtree_3.14.0                    
[105] bit64_4.6.0-1                     RcppAnnoy_0.0.22                 
[107] GenomeInfoDb_1.42.3               bslib_0.9.0                      
[109] irlba_2.3.5.1                     KernSmooth_2.23-20               
[111] seqLogo_1.72.0                    SeuratDisk_0.0.0.9021            
[113] colorspace_2.1-1                  BiocGenerics_0.52.0              
[115] DBI_1.2.3                         tidyselect_1.2.1                 
[117] bit_4.6.0                         compiler_4.4.2                   
[119] curl_6.4.0                        hdf5r_1.3.12                     
[121] DelayedArray_0.32.0               plotly_4.11.0                    
[123] rtracklayer_1.66.0                caTools_1.18.3                   
[125] lmtest_0.9-40                     rappdirs_0.3.3                   
[127] digest_0.6.37                     goftest_1.2-3                    
[129] presto_1.0.0                      spatstat.utils_3.1-4             
[131] rmarkdown_2.29                    RhpcBLASctl_0.23-42              
[133] XVector_0.46.0                    htmltools_0.5.8.1                
[135] pkgconfig_2.0.3                   MatrixGenerics_1.18.1            
[137] fastmap_1.2.0                     ensembldb_2.30.0                 
[139] rlang_1.1.6                       htmlwidgets_1.6.4                
[141] UCSC.utils_1.2.0                  shiny_1.11.1                     
[143] farver_2.1.2                      jquerylib_0.1.4                  
[145] zoo_1.8-14                        jsonlite_2.0.0                   
[147] BiocParallel_1.40.2               GOSemSim_2.32.0                  
[149] R.oo_1.27.1                       RCurl_1.98-1.17                  
[151] magrittr_2.0.3                    GenomeInfoDbData_1.2.13          
[153] ggplotify_0.1.2                   dotCall64_1.2                    
[155] Rhdf5lib_1.28.0                   ape_5.8-1                        
[157] babelgene_22.9                    reticulate_1.42.0                
[159] stringi_1.8.7                     zlibbioc_1.52.0                  
[161] MASS_7.3-58.2                     plyr_1.8.9                       
[163] parallel_4.4.2                    listenv_0.9.1                    
[165] ggrepel_0.9.6                     CNEr_1.42.0                      
[167] deldir_2.0-4                      Biostrings_2.74.1                
[169] splines_4.4.2                     tensor_1.5.1                     
[171] hms_1.1.3                         BSgenome.Hsapiens.UCSC.hg38_1.4.5
[173] igraph_2.1.4                      spatstat.geom_3.4-1              
[175] RcppHNSW_0.6.0                    reshape2_1.4.4                   
[177] stats4_4.4.2                      TFMPvalue_0.0.9                  
[179] XML_3.99-0.18                     evaluate_1.0.4                   
[181] tzdb_0.5.0                        JASPAR2020_0.99.10               
[183] httpuv_1.6.16                     RANN_2.6.2                       
[185] purrr_1.0.4                       polyclip_1.10-7                  
[187] SeuratData_0.2.2.9002             scattermore_1.2                  
[189] xtable_1.8-4                      restfulr_0.0.16                  
[191] AnnotationFilter_1.30.0           RSpectra_0.16-2                  
[193] tidytree_0.4.6                    later_1.4.2                      
[195] googledrive_2.1.1                 viridisLite_0.4.2                
[197] tibble_3.3.0                      aplot_0.2.8                      
[199] memoise_2.0.1                     AnnotationDbi_1.68.0             
[201] GenomicAlignments_1.42.0          IRanges_2.40.1                   
[203] cluster_2.1.4                     globals_0.18.0                   
LS0tCnRpdGxlOiAiSW50ZWdyYXRlIFZpcnVzIHNjUk5BLXNlcSBkYXRhIgojIDA2LTI3LTIwMjUKIyB1c2VkIHNjUk5BLXNlcSBkYXRhOiBFX01UQUJfMTIzOTIgYW5kIEdTRTE2NDM4MQpvdXRwdXQ6IGh0bWxfbm90ZWJvb2sKLS0tCgoKYGBge3J9CmxpYnJhcnkoU2V1cmF0KQpsaWJyYXJ5KGdncGxvdDIpCmxpYnJhcnkoZHBseXIpCmxpYnJhcnkoaGFybW9ueSkKbGlicmFyeShjbHVzdGVyUHJvZmlsZXIpCmxpYnJhcnkoZW5yaWNocGxvdCkKbGlicmFyeShtc2lnZGJyKQpsaWJyYXJ5KHN0cmluZ3IpCmxpYnJhcnkodGlkeXIpCmxpYnJhcnkocGF0Y2h3b3JrKQpsaWJyYXJ5KEF6aW11dGgpCnNldC5zZWVkKDIwMjUpCmBgYAoKIyBsb2FkIGRhdGEKYGBge3J9Cm91dF9wYXRoIDwtICJ+L3Byb2plY3RzLzIwMjQvUkEvRmx1X2ludGVncmF0aW9uLyIKIyBsb2FkIEVfTVRBQl8xMjM5Mgppbl9wYXRoIDwtICJ+L3Byb2plY3RzLzIwMjQvUkEvRV9NVEFCXzEyMzkyL291dC8iCkVfTVRBQl8xMjM5Ml9zY2RhdGEgPC0gcmVhZFJEUyhwYXN0ZTAoaW5fcGF0aCwiRV9NVEFCXzEyMzkyX3NjZGF0YS5yZHMiKSkKRV9NVEFCXzEyMzkyX21ldGFkYXRhIDwtIHJlYWRSRFMocGFzdGUwKGluX3BhdGgsIkVfTVRBQl8xMjM5Ml9tZXRhZGF0YS5yZHMiKSkKbGVuZ3RoKEVfTVRBQl8xMjM5Ml9zY2RhdGEpCiMgbG9hZCBHU0UxNjQzODEKaW5fcGF0aCA8LSAifi9wcm9qZWN0cy8yMDI0L1JBL0dTRTE2NDM4MV9SQVcvb3V0LyIKR1NFMTY0MzgxX3NjZGF0YSA8LSByZWFkUkRTKHBhc3RlMChpbl9wYXRoLCJHU0UxNjQzODFfc2NkYXRhLnJkcyIpKQpHU0UxNjQzODFfbWV0YWRhdGEgPC0gcmVhZFJEUyhwYXN0ZTAoaW5fcGF0aCwiR1NFMTY0MzgxX21ldGFkYXRhLnJkcyIpKQpsZW5ndGgoR1NFMTY0MzgxX3NjZGF0YSkKYGBgCgojIHNlbGVjdCBwYXNzZWQgY3V0b2ZmIGNlbGxzIGluIGVhY2ggY29ob3J0CmBgYHtyfQojIHNhbXBsZSBuYW1lCmRfbmFtZV8xIDwtIHVuaXF1ZShFX01UQUJfMTIzOTJfbWV0YWRhdGEkb3JpZy5pZGVudCkKZF9uYW1lXzIgPC0gdW5pcXVlKEdTRTE2NDM4MV9tZXRhZGF0YSRvcmlnLmlkZW50KQpkX25hbWUgPC0gYyhkX25hbWVfMSxkX25hbWVfMikKIyBjZWxsIGxpc3QKY2VsbF9saXN0XzEgPC0gcm93bmFtZXMoRV9NVEFCXzEyMzkyX21ldGFkYXRhKQpjZWxsX2xpc3RfMiA8LSByb3duYW1lcyhHU0UxNjQzODFfbWV0YWRhdGEpCmhlYWQoY2VsbF9saXN0XzEpCmhlYWQoY2VsbF9saXN0XzIpCmNlbGxfbGlzdCA8LSBjKGNlbGxfbGlzdF8xLCBjZWxsX2xpc3RfMikKIyBtZXJnZQpzY2RhdGEgPC0gYyhFX01UQUJfMTIzOTJfc2NkYXRhLCBHU0UxNjQzODFfc2NkYXRhKQpTb2JqIDwtIG1lcmdlKHNjZGF0YVtbMV1dLCB5ID0gc2NkYXRhWy0xXSwgYWRkLmNlbGwuaWRzID0gZF9uYW1lKQpTb2JqW1sicGVyY2VudC5tdCJdXSA8LSBQZXJjZW50YWdlRmVhdHVyZVNldChTb2JqLCBwYXR0ZXJuID0gIl5NVC0iKQojIHNlbGVjdCBwYXNzZWQgY2VsbHMgClNvYmogPC0gc3Vic2V0KFNvYmosIGNlbGxzID0gY2VsbF9saXN0KQpsZW5ndGgoY2VsbF9saXN0KQpkaW0oU29iaikKYGBgCgojIG5vcm1hbGl6ZSwgc2NhbGUsIFBDQQpgYGB7cn0KU29iaiA8LSBOb3JtYWxpemVEYXRhKFNvYmosIG5vcm1hbGl6YXRpb24ubWV0aG9kID0gIkxvZ05vcm1hbGl6ZSIsIHNjYWxlLmZhY3RvciA9IDEwMDAwKQpTb2JqIDwtIEZpbmRWYXJpYWJsZUZlYXR1cmVzKFNvYmosIHNlbGVjdGlvbi5tZXRob2QgPSAidnN0IiwgbmZhdHVyZXMgPSAyMDAwKQpTb2JqIDwtIFNjYWxlRGF0YShTb2JqKQpTb2JqIDwtIFJ1blBDQShTb2JqKQpgYGAKCiMgc2VsZWN0IGN1dG9mZiBvZiBkaW1lbnNpb24uIApgYGB7cn0KRWxib3dQbG90KFNvYmopCmBgYAojIFBDQSBzcGFjZQpgYGB7cn0KRGltUGxvdChTb2JqLCBncm91cC5ieSA9ICJvcmlnLmlkZW50IikKYGBgCiMgVU1BUApgYGB7cn0KU29iaiA8LSBSdW5VTUFQKFNvYmosIGRpbXMgPSAxOjE1LCByZWR1Y3Rpb24gPSAicGNhIiwgcmVkdWN0aW9uLm5hbWUgPSAidW1hcC51bmludGVncmF0ZWQiKQpwMCA8LSBEaW1QbG90KFNvYmosIHJlZHVjdGlvbiA9ICJ1bWFwLnVuaW50ZWdyYXRlZCIsIGdyb3VwLmJ5ID0gYygib3JpZy5pZGVudCIsICJCYXRjaCIpKQpwMApgYGAKCiMgSGFybW9ueSBiYXRjaCBjb3JyZWN0aW9uCmBgYHtyfQpTb2JqIDwtIFJ1bkhhcm1vbnkoU29iaiwib3JpZy5pZGVudCIsIHBsb3RfY29udmVyZ2VuY2UgPSBULCB0aGV0YSA9IDMpCmhhcm1vbnlfZW1iZWRkaW5ncyA8LSBFbWJlZGRpbmdzKFNvYmosICJoYXJtb255IikKaGFybW9ueV9lbWJlZGRpbmdzWzE6NSwgMTo1XQpEaW1QbG90KFNvYmosIHJlZHVjdGlvbiA9ICJoYXJtb255IiwgZ3JvdXAuYnkgPSAib3JpZy5pZGVudCIpCmBgYAojIFVNQVAgd2l0aCBoYXJtb255CmBgYHtyfQpTb2JqIDwtIFJ1blVNQVAoU29iaiwgZGltcyA9IDE6MTUsIHJlZHVjdGlvbiA9ICJoYXJtb255IikKYGBgCgpgYGB7cn0KRGltUGxvdChTb2JqLCByZWR1Y3Rpb24gPSAidW1hcCIsIGdyb3VwLmJ5ID0gIm9yaWcuaWRlbnQiKQpgYGAKIyBDaGVjayB3aGV0aGVyIHRoZSBoYXJtb255IGJhdGNoIGNvcnJlY3Rpb24gd29yayB3ZWxsIG9yIG5vdC4gCmBgYHtyfQpJZGVudHMoU29iaikgPC0gU29iaiRvcmlnLmlkZW50CnJlbmFtZV9tYXAgPC0gYygKICAgIHNldE5hbWVzKHJlcCgiRV9NVEFCXzEyMzkyIiwgbGVuZ3RoKGRfbmFtZV8xKSksIGRfbmFtZV8xKSwKICAgIHNldE5hbWVzKHJlcCgiR1NFMTY0MzgxIiwgbGVuZ3RoKGRfbmFtZV8yKSksIGRfbmFtZV8yKQopClNvYmogPC0gUmVuYW1lSWRlbnRzKFNvYmosIHJlbmFtZV9tYXApCnVuaXF1ZShJZGVudHMoU29iaikpClNvYmokYmF0Y2ggPC0gSWRlbnRzKFNvYmopCmBgYAoKYGBge3J9CkRpbVBsb3QoU29iaiwgZ3JvdXAuYnkgPSAiYmF0Y2giKQpEaW1QbG90KFNvYmosIHJlZHVjdGlvbiA9ICJ1bWFwIiwgZ3JvdXAuYnkgPSAiYmF0Y2giKQoKYGBgCgpgYGB7cn0KU29iagpTb2JqIDwtIEZpbmROZWlnaGJvcnMoU29iaiwgZGltcyA9IDE6MTUsIHJlZHVjdGlvbiA9ICJoYXJtb255IikKU29iaiA8LSBGaW5kQ2x1c3RlcnMoU29iaiwgcmVkdWN0aW9uID0gInVtYXAiLCByZXNvbHV0aW9uID0gMSkKU29iakByZWR1Y3Rpb25zCmBgYAoKYGBge3J9CkRpbVBsb3QoU29iaiwgcmVkdWN0aW9uID0gInVtYXAiLCBncm91cC5ieSA9ICJSTkFfc25uX3Jlcy4xIikKYGBgCgpgYGB7cn0KU29iaiA8LSBKb2luTGF5ZXJzKFNvYmopCmBgYAoKYGBge3J9CkRpbVBsb3QoU29iaiwgcmVkdWN0aW9uID0gInVtYXAiLCBsYWJlbCA9IFQpCmBgYAoKIyBsb2FkIHJlZmVyZW5jZSBkYXRhCmBgYHtyfQpTTEVfU29iaiA8LSByZWFkUkRTKCJ+L3Byb2plY3RzLzIwMjQvUkEvSW50ZWdyYXRpb25fU0xFX1NvYmoucmRzIikKU0xFX1NvYmogPC0gUnVuVU1BUChTTEVfU29iaiwgcmVkdWN0aW9uID0gImhhcm1vbnkiLGRpbXMgPSAxOjE1LCByZXR1cm4ubW9kZWwgPSBUUlVFKQpgYGAKCmBgYHtyfQpEaW1QbG90KFNMRV9Tb2JqLCByZWR1Y3Rpb24gPSAidW1hcCIsIGxhYmVsID0gVCwgcmVwZWwgPSBUKQpEaW1QbG90KFNMRV9Tb2JqLCByZWR1Y3Rpb24gPSAidW1hcC51bmludGVncmF0ZWQiKQoKYGBgCgoKYGBge3J9CiMgaWYgKCFmaWxlLmV4aXN0cyhwYXN0ZTAob3V0X3BhdGgsIkludGVncmF0aW9uX1ZpcnVzLlJEYXRhIikpKSB7CiMgICBzYXZlLmltYWdlKGZpbGUgPSBwYXN0ZTAob3V0X3BhdGgsIkludGVncmF0aW9uX1ZpcnVzLlJEYXRhIikpCiMgfQpgYGAKCgojIGNhbGN1bGF0ZSB0cmFuc2ZlciBhbmNob3JzCmBgYHtyfQphbmNob3JzIDwtIEZpbmRUcmFuc2ZlckFuY2hvcnMocmVmZXJlbmNlID0gU0xFX1NvYmosCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBxdWVyeSA9IFNvYmosCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICBub3JtYWxpemF0aW9uLm1ldGhvZCA9ICJMb2dOb3JtYWxpemUiLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgcmVmZXJlbmNlLnJlZHVjdGlvbiA9ICJoYXJtb255IiwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIGRpbXMgPSAxOjE1KQpgYGAKCiMgTWFwcGluZwpgYGB7cn0KU29iaiA8LSBNYXBRdWVyeShhbmNob3JzZXQgPSBhbmNob3JzLAogICAgICAgICAgICAgICAgIHF1ZXJ5ID0gU29iaiwKICAgICAgICAgICAgICAgICByZWZlcmVuY2UgPSBTTEVfU29iaiwKICAgICAgICAgICAgICAgICByZWZkYXRhID0gbGlzdChjZWxsdHlwZSA9ICJhbm5vX3RleHQiKSwKICAgICAgICAgICAgICAgICByZWZlcmVuY2UucmVkdWN0aW9uID0gImhhcm1vbnkiLAogICAgICAgICAgICAgICAgIHJlZHVjdGlvbi5tb2RlbCA9ICJ1bWFwIikKYGBgCgpgYGB7cn0KSWRlbnRzKFNvYmopIDwtIFNvYmokcHJlZGljdGVkLmNlbGx0eXBlCklkZW50cyhTb2JqKSA8LSBmYWN0b3IoeCA9IElkZW50cyhTb2JqKSwgbGV2ZWxzID0gYygiTmFpdmUgQiIsICJBY3RpdmF0ZWQgQiIsICJVbnN3aXRjaGVkIG1lbW9yeSBCIiwgIkNEMjcrIG1lbW9yeSBCIiwgIlBsYXNtYWJsYXN0IiwgIlBsYXNtYSIsICJUIGNlbGwgbGlrZSBCIiwgIk15ZWxvaWQgbGlrZSBCIiwgIk1lZ2FrYXJ5b2N5dGUgbGlrZSBCIikpCmxldmVscyhJZGVudHMoU29iaikpCmxldmVscyhTTEVfU29iaiRhbm5vX3RleHQpClNvYmokcHJlZGljdGVkLmNlbGx0eXBlIDwtIElkZW50cyhTb2JqKQpgYGAKCmBgYHtyfQp0YWJsZShTb2JqJHByZWRpY3RlZC5jZWxsdHlwZSkKYGBgCgpgYGB7cn0KSWRlbnRzKFNvYmopIDwtIFNvYmokcHJlZGljdGVkLmNlbGx0eXBlClNvYmogPC0gc3Vic2V0KFNvYmosIGlkZW50cyA9IHNldGRpZmYobGV2ZWxzKFNvYmopLCAiTWVnYWthcnlvY3l0ZSBsaWtlIEIiKSkKdGFibGUoU29iaiRwcmVkaWN0ZWQuY2VsbHR5cGUpCmBgYAoKCmBgYHtyfQpJZGVudHMoU29iaikgPC0gU29iaiRSTkFfc25uX3Jlcy4xCiMgcDEgPC0gRGltUGxvdChTb2JqLCByZWR1Y3Rpb24gPSAidW1hcCIsIGdyb3VwLmJ5ID0gIlJOQV9zbm5fcmVzLjEiKQpwMiA8LSBEaW1QbG90KFNvYmosIHJlZHVjdGlvbiA9ICJ1bWFwIiwgbGFiZWwgPSBUKQpJZGVudHMoU29iaikgPC0gU29iaiRwcmVkaWN0ZWQuY2VsbHR5cGUKcDMgPC0gRGltUGxvdChTb2JqLCBsYWJlbCA9IFRSVUUsIHJlcGVsID0gVFJVRSwgcmVkdWN0aW9uID0gInJlZi51bWFwIikgKyBOb0xlZ2VuZCgpCnA0IDwtIERpbVBsb3QoU29iaiwgbGFiZWwgPSBUUlVFLCByZXBlbCA9IFRSVUUsIHJlZHVjdGlvbiA9ICJyZWYudW1hcCIpCnA1IDwtIERpbVBsb3QoU29iaiwgbGFiZWwgPSBUUlVFLCByZXBlbCA9IFRSVUUsIHJlZHVjdGlvbiA9ICJ1bWFwIikgKyBOb0xlZ2VuZCgpCnA2IDwtIERpbVBsb3QoU29iaiwgbGFiZWwgPSBUUlVFLCByZXBlbCA9IFRSVUUsIHJlZHVjdGlvbiA9ICJ1bWFwIikKYGBgCgoKCgpgYGB7cn0KRV9NVEFCXzEyMzkyX21ldGFkYXRhJGFubm9fQUJDX0xSQiA8LSBOQQoKRV9NVEFCXzEyMzkyX21ldGFkYXRhJGFubm9fQUJDX0xSQltFX01UQUJfMTIzOTJfbWV0YWRhdGEkTFNQMV9zY29yZV9wb3M9PVQgJiBFX01UQUJfMTIzOTJfbWV0YWRhdGEkQUJDX3Njb3JlX3Bvcz09Rl0gPC0gIkxSQiIKRV9NVEFCXzEyMzkyX21ldGFkYXRhJGFubm9fQUJDX0xSQltFX01UQUJfMTIzOTJfbWV0YWRhdGEkTFNQMV9zY29yZV9wb3M9PUYgJiBFX01UQUJfMTIzOTJfbWV0YWRhdGEkQUJDX3Njb3JlX3Bvcz09VF0gPC0gIkFCQyIKRV9NVEFCXzEyMzkyX21ldGFkYXRhJGFubm9fQUJDX0xSQltFX01UQUJfMTIzOTJfbWV0YWRhdGEkTFNQMV9zY29yZV9wb3M9PVQgJiBFX01UQUJfMTIzOTJfbWV0YWRhdGEkQUJDX3Njb3JlX3Bvcz09VF0gPC0gIkxSQiAmIEFCQyIKRV9NVEFCXzEyMzkyX21ldGFkYXRhJGFubm9fQUJDX0xSQltFX01UQUJfMTIzOTJfbWV0YWRhdGEkTFNQMV9zY29yZV9wb3M9PUYgJiBFX01UQUJfMTIzOTJfbWV0YWRhdGEkQUJDX3Njb3JlX3Bvcz09Rl0gPC0gIkNsYXNzaWNhbCBCIgoKdW5pcXVlKEVfTVRBQl8xMjM5Ml9tZXRhZGF0YSRhbm5vX0FCQ19MUkIpCnVuaXF1ZShHU0UxNjQzODFfbWV0YWRhdGEkYW5ub19BQkNfTFJCKQoKYW5ub192YXJzIDwtIGMoImFubm9fQUJDX0xSQiIsICJMU1AxX3Njb3JlIiwgIkFCQ19zY29yZSIsICJMU1AxX3Njb3JlX3BvcyIsICJBQkNfc2NvcmVfcG9zIikKCmZvciAodmFyIGluIGFubm9fdmFycykgewogIGwxIDwtIEVfTVRBQl8xMjM5Ml9tZXRhZGF0YVtbdmFyXV0KICBsMiA8LSBHU0UxNjQzODFfbWV0YWRhdGFbW3Zhcl1dCiAgCiAgaWYgKGlzLmZhY3RvcihsMSkpIGwxIDwtIGFzLmNoYXJhY3RlcihsMSkKICBpZiAoaXMuZmFjdG9yKGwyKSkgbDIgPC0gYXMuY2hhcmFjdGVyKGwyKQogIAogIG5hbWVzKGwxKSA8LSByb3duYW1lcyhFX01UQUJfMTIzOTJfbWV0YWRhdGEpCiAgbmFtZXMobDIpIDwtIHJvd25hbWVzKEdTRTE2NDM4MV9tZXRhZGF0YSkKICAKICBjb21iaW5lZF9hbm5vIDwtIGMobDEsIGwyKQogIAogIFNvYmpbW3Zhcl1dIDwtIGNvbWJpbmVkX2Fubm9bY29sbmFtZXMoU29iaildCn0Kcm0obDEsIGwyLCBjb21iaW5lZF9hbm5vKQpgYGAKCgpgYGB7cn0KbXlfY29sb3JzIDwtIGMoIkxSQiIgPSAiI2ZmMDAwMCIsICJBQkMiID0gIiNmZmZmMDAiLCAiTFJCICYgQUJDIiA9ICIjMDBiMGZmIiwgIkNsYXNzaWNhbCBCIiA9ICIjQkZCRkJGIikKSWRlbnRzKFNvYmopIDwtIFNvYmokYW5ub19BQkNfTFJCCnA3IDwtIERpbVBsb3QoU29iaiwgcmVkdWN0aW9uID0gInVtYXAiLCBjb2xzID0gbXlfY29sb3JzKSArIE5vTGVnZW5kKCkKRGltUGxvdChTb2JqLCByZWR1Y3Rpb24gPSAidW1hcCIsIGNvbHMgPSBteV9jb2xvcnMpCmBgYAoKYGBge3J9CmNlbGx0eXBlLnByb3AgPC0gYXMuZGF0YS5mcmFtZSh0YWJsZShTb2JqJHByZWRpY3RlZC5jZWxsdHlwZSwgU29iaiRhbm5vX0FCQ19MUkIpKQojIGhlYWQoY2VsbHR5cGUucHJvcCkKY2VsbHR5cGUucHJvcCA8LSBjZWxsdHlwZS5wcm9wICU+JQogICAgcmVuYW1lKGNlbGx0eXBlID0gVmFyMSwgYEFCQyBMUkJgID0gVmFyMikgJT4lCiAgICBncm91cF9ieShjZWxsdHlwZSkgJT4lCiAgICBtdXRhdGUoUGVyY2VudCA9IEZyZXEgLyBzdW0oRnJlcSkqMTAwKQpnZ3Bsb3QoY2VsbHR5cGUucHJvcCwgYWVzKHggPSBjZWxsdHlwZSwgeSA9IFBlcmNlbnQsIGZpbGwgPSBgQUJDIExSQmApKSsKICAgIGdlb21fYmFyKHN0YXQgPSAiaWRlbnRpdHkiLCB3aWR0aCA9IDAuOSkgKyAKICAgIHNjYWxlX2ZpbGxfbWFudWFsKHZhbHVlcyA9IG15X2NvbG9ycykgKwogICAgdGhlbWUoYXhpcy50ZXh0LnggPSBlbGVtZW50X3RleHQoYW5nbGUgPSA0NSwgdmp1c3QgPSAxLCBoanVzdCA9IDEpLCBheGlzLnRpdGxlLnggPSBlbGVtZW50X2JsYW5rKCksIHBsb3QubWFyZ2luID0gbWFyZ2luKDEwLCAxMCwgMjAsIDEwKSkKCgpjZWxsdHlwZV9wZXJjZW50X3dpZGUgPC0gY2VsbHR5cGUucHJvcCAlPiUKICBzZWxlY3QoY2VsbHR5cGUsIGBBQkMgTFJCYCwgUGVyY2VudCkgJT4lCiAgcGl2b3Rfd2lkZXIobmFtZXNfZnJvbSA9IGBBQkMgTFJCYCwgdmFsdWVzX2Zyb20gPSBQZXJjZW50KQpoZWFkKGNlbGx0eXBlX3BlcmNlbnRfd2lkZSkKCmNlbGx0eXBlX2ZyZXFfd2lkZSA8LSBjZWxsdHlwZS5wcm9wICU+JQogIHNlbGVjdChjZWxsdHlwZSwgYEFCQyBMUkJgLCBGcmVxKSAlPiUKICBwaXZvdF93aWRlcihuYW1lc19mcm9tID0gYEFCQyBMUkJgLCB2YWx1ZXNfZnJvbSA9IEZyZXEpCmhlYWQoY2VsbHR5cGVfZnJlcV93aWRlKQp3cml0ZS5jc3YoY2VsbHR5cGVfcGVyY2VudF93aWRlLCBwYXN0ZTAob3V0X3BhdGgsIm91dC92aXJ1c19iYXJwbG90X3BlcmNlbnQuY3N2IikpCndyaXRlLmNzdihjZWxsdHlwZV9mcmVxX3dpZGUsIHBhc3RlMChvdXRfcGF0aCwib3V0L3ZpcnVzX2JhcnBsb3RfZnJlcS5jc3YiKSkKYGBgCgoKYGBge3J9CnAyCnAzCnA0CnA1CnA2CnA3CnBkZihwYXN0ZTAob3V0X3BhdGgsICJBemltdXRoX3ByZWRpY3RlZF91bWFwXzA2MzAyNS5wZGYiKSwgd2lkdGggPSA0LCBoZWlnaHQgPSA0KQpwMgpwMwpwNApwNQpwNgpwNwpkZXYub2ZmKCkKYGBgCgpgYGB7cn0KSWRlbnRzKFNvYmopIDwtIFNvYmokcHJlZGljdGVkLmNlbGx0eXBlCnA2IDwtIERpbVBsb3QoU29iaiwgbGFiZWwgPSBUUlVFLCByZXBlbCA9IFRSVUUsIHJlZHVjdGlvbiA9ICJ1bWFwIikKcGRmKHBhc3RlMChvdXRfcGF0aCwgIm5ld192aXJ1c191bWFwXzA4MDQyNS5wZGYiKSwgd2lkdGggPSA0LCBoZWlnaHQgPSA0KQpwNgpwNwpkZXYub2ZmKCkKcDZfMSA8LSBEaW1QbG90KFNvYmosIGxhYmVsID0gRkFMU0UsIHJlcGVsID0gVFJVRSwgcmVkdWN0aW9uID0gInVtYXAiKSArCiAgbGFicyh4ID0gIiIsIHkgPSAiIikgKwogIHRoZW1lKAogICAgbGVnZW5kLnBvc2l0aW9uID0gIm5vbmUiLAogICAgYXhpcy50ZXh0LnggPSBlbGVtZW50X2JsYW5rKCksCiAgICBheGlzLnRleHQueSA9IGVsZW1lbnRfYmxhbmsoKSwKICAgIGF4aXMudGlja3MueCA9IGVsZW1lbnRfbGluZSgpLAogICAgYXhpcy50aWNrcy55ID0gZWxlbWVudF9saW5lKCkKICAgICkKdGlmZihmaWxlbmFtZSA9IHBhc3RlMChvdXRfcGF0aCwidW1hcF92aXJ1cy50aWYiKSwgd2lkdGggPSA0LCBoZWlnaHQgPSA0LCB1bml0cyA9ICJpbiIsIHJlcyA9IDMwMCwgY29tcHJlc3Npb24gPSAibHp3IikKcHJpbnQocDZfMSkKZGV2Lm9mZigpCgpwN18xIDwtIHA3ICsKICBsYWJzKHggPSAiIiwgeSA9ICIiKSArCiAgdGhlbWUoCiAgICBsZWdlbmQucG9zaXRpb24gPSAibm9uZSIsCiAgICBheGlzLnRleHQueCA9IGVsZW1lbnRfYmxhbmsoKSwKICAgIGF4aXMudGV4dC55ID0gZWxlbWVudF9ibGFuaygpLAogICAgYXhpcy50aWNrcy54ID0gZWxlbWVudF9saW5lKCksCiAgICBheGlzLnRpY2tzLnkgPSBlbGVtZW50X2xpbmUoKQogICAgKQp0aWZmKGZpbGVuYW1lID0gcGFzdGUwKG91dF9wYXRoLCJ1bWFwX0xSQl9BQkMudGlmIiksIHdpZHRoID0gNCwgaGVpZ2h0ID0gNCwgdW5pdHMgPSAiaW4iLCByZXMgPSAzMDAsIGNvbXByZXNzaW9uID0gImx6dyIpCnByaW50KHA3XzEpCmRldi5vZmYoKQoKYGBgCgoKYGBge3J9CnRhYmxlKFNvYmokYW5ub19BQkNfTFJCKQpgYGAKCmBgYHtyfQpsaWJyYXJ5KGdncG9pbnRkZW5zaXR5KQpsaWJyYXJ5KHNjYWxlcykKYGBgCgpgYGB7cn0KZGYgPC0gRmV0Y2hEYXRhKFNvYmosIHZhcnMgPSBjKCJ1bWFwXzEiLCAidW1hcF8yIiwgIkxTUDFfc2NvcmVfcG9zIiwgIkFCQ19zY29yZV9wb3MiLCAiYW5ub19BQkNfTFJCIikpCmRmX0xTUDEgPC0gZGYgJT4lIGZpbHRlcihMU1AxX3Njb3JlX3BvcykKZGZfQUJDIDwtIGRmICU+JSBmaWx0ZXIoQUJDX3Njb3JlX3BvcykKYGBgCgpgYGB7cn0KdDEgPC0gZ2dwbG90KGRmX0xTUDEsIGFlcyh4ID0gdW1hcF8xLCB5ID0gdW1hcF8yKSkgKwogIGdlb21fcG9pbnRkZW5zaXR5KG1ldGhvZCA9ICJuZWlnaGJvcnMiLCBzaXplID0gMC4xKSArCiAgc2NhbGVfY29sb3JfZ3JhZGllbnRuKAogICAgY29sb3JzID0gYygiZ3JleSIsICJncmV5IiwgInJlZCIpLAogICAgdmFsdWVzID0gcmVzY2FsZShjKDAsIDM1MCwgMTIwMCkpLAogICAgbGltaXRzID0gYygwLCAxMDAwKQogICkgKwogICAgdGhlbWVfY2xhc3NpYygpCgp0MiA8LSBnZ3Bsb3QoZGZfQUJDLCBhZXMoeCA9IHVtYXBfMSwgeSA9IHVtYXBfMikpICsKICBnZW9tX3BvaW50ZGVuc2l0eShtZXRob2QgPSAibmVpZ2hib3JzIiwgc2l6ZSA9IDAuMSkgKwogIHNjYWxlX2NvbG9yX2dyYWRpZW50bigKICAgIGNvbG9ycyA9IGMoImdyZXkiLCAiZ3JleSIsICJyZWQiKSwKICAgIHZhbHVlcyA9IHJlc2NhbGUoYygwLCAzNTAsIDEyMDApKSwKICAgIGxpbWl0cyA9IGMoMCwgMTAwMCkKICApICsKICAgIHRoZW1lX2NsYXNzaWMoKQp0MQp0MgpgYGAKCgoKCmBgYHtyfQpwZGYocGFzdGUwKG91dF9wYXRoLCJMUkJfQUJDX3VtYXBfZGVuc2l0eV8wODA0MjUucGRmIiksIHdpZHRoID0gNiwgaGVpZ2h0ID0gNSkKdDEKdDIKZGV2Lm9mZigpCmBgYAoKYGBge3J9CnQxXzEgPC0gdDEgKwpsYWJzKHggPSAiIiwgeSA9ICIiKSArCnRoZW1lKAogIGxlZ2VuZC5wb3NpdGlvbiA9ICJub25lIiwKICBheGlzLnRleHQueCA9IGVsZW1lbnRfYmxhbmsoKSwKICBheGlzLnRleHQueSA9IGVsZW1lbnRfYmxhbmsoKSwKICBheGlzLnRpY2tzLnggPSBlbGVtZW50X2xpbmUoKSwKICBheGlzLnRpY2tzLnkgPSBlbGVtZW50X2xpbmUoKQogICkKdDJfMSA8LSB0MiArCmxhYnMoeCA9ICIiLCB5ID0gIiIpICsKdGhlbWUoCiAgbGVnZW5kLnBvc2l0aW9uID0gIm5vbmUiLAogIGF4aXMudGV4dC54ID0gZWxlbWVudF9ibGFuaygpLAogIGF4aXMudGV4dC55ID0gZWxlbWVudF9ibGFuaygpLAogIGF4aXMudGlja3MueCA9IGVsZW1lbnRfbGluZSgpLAogIGF4aXMudGlja3MueSA9IGVsZW1lbnRfbGluZSgpCiAgKQp0MV8xCnQyXzEKYGBgCgpgYGB7cn0KdGlmZihmaWxlbmFtZSA9IHBhc3RlMChvdXRfcGF0aCwiTFJCX2RlbnNpdHlfdW1hcC50aWYiKSwgd2lkdGggPSA0LCBoZWlnaHQgPSA0LCB1bml0cyA9ICJpbiIsIHJlcyA9IDMwMCwgY29tcHJlc3Npb24gPSAibHp3IikKcHJpbnQodDFfMSkKZGV2Lm9mZigpCgp0aWZmKGZpbGVuYW1lID0gcGFzdGUwKG91dF9wYXRoLCJBQkNfZGVuc2l0eV91bWFwLnRpZiIpLCB3aWR0aCA9IDQsIGhlaWdodCA9IDQsIHVuaXRzID0gImluIiwgcmVzID0gMzAwLCBjb21wcmVzc2lvbiA9ICJsenciKQpwcmludCh0Ml8xKQpkZXYub2ZmKCkKYGBgCgpgYGB7cn0KaWYgKCFmaWxlLmV4aXN0cyhwYXN0ZTAob3V0X3BhdGgsInZpcnVzX2ludGVncmF0ZWRfZGF0YV8wODA0MjUucmRzIikpKSB7CiAgc2F2ZVJEUyhTb2JqLCBmaWxlID0gcGFzdGUwKG91dF9wYXRoLCJ2aXJ1c19pbnRlZ3JhdGVkX2RhdGFfMDgwNDI1LnJkcyIpKQp9CmBgYAoKCgpgYGB7cn0Kc2Vzc2lvbkluZm8oKQpgYGAKCgoKCgoKCgoKCgoKCg==